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\  ABSTRACT 

one-dimensional  numerical  code  has  been  written,  which  allows  the 
fj  A  time-dependent  simulation  of  the  Child-Langmuir  diode  and  the  plasma  sheath 
/  proBftmr-Our  model  assumes  a  semi-infinite  plasma  reservoir  from  which  parti¬ 
cles  can  be  emitted  into  a  vacuum  region  where  they  can  be  accelerated  or 
decelerated  by  means  of  an  external  potential  on  the  anode,  until  they  are 
eventually  absorbed  by  the  anode  or  sent  back  to  the  reservoir. 

Physical  results  are  shown  for  both  the  time-independent  and  the  time- 
dependent  cases.  Numerical  methods  are  presented  for  the  evaluation  of 
relevant  physical  quantities.  The  numerical  effects  are  discussed.  ^ 
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1.  INTRODUCTION 

The  growing  interest  in  many  different  disciplines  for  diode-like  problems,  along  with  the 
highly  nonlinear  behavior  of  such  configurations,  have  stimulated  the  search  for  solutions  to 
these  problems  by  means  of  numerical  simulation. 

Scientists  in  the  fields  of  plasma  physics  and  space  sciences  are  especially  interested  in  the 
results  of  such  simulations,  since  our  problem  has  direct  applications  in  the  domains  of  Q- 
machines,  plasma-surface  interactions  and  effects  of  the  first  wall  on  the  plasma  within  a  fusion 
reactor,  plasma  insulators,  Langmuir  probes,  probe-wall  effects,  etc. 

The  general  purpose  of  the  study  is  to  simulate  the  behavior  of  the  plasma  sheath  by 
means  of  a  unidimensional  numerical  code,  in  order  to  reach  by  simulation  the  essence  of  the 
physics  of  this  problem.  Our  goal  in  this  project  was  to  build  and  test  the  tool  (a  numerical 
code)  necessary  for  such  a  search. 

Many  people  have  studied  diode-like  problems  analytically,  experimentally  and  computa¬ 
tionally:  see  Child  (1911),  Langmuir  (1923),  Pierce  (1948),  Tien  St  Moshman  (1956),  Self 
(1963),  Birdsall  St  Bridges  (1966),  Kuhn  (1979).  In  their  classic  paper,  I.  Langmuir  and  H. 
Mott-Smith  (1924)  studied  the  time-independent  solution  of  the  diode  problem.  However,  the 
stability  of  this  solution  was  not  described  and  several  questions  had  aot  been  answered:  was 
there  some  kind  of  oscillating  mode,  or  was  the  equilibrium  perfectly  time-independent?  And 
in  case  such  an  oscillation  existed,  was  the  time-independent  solution  just  the  time-average  of 
the  oscillation? 

In  our  study,  we  have  been  especially  interested  in  measuring  the  numerical  effects  due  to 
the  model  we  chose;  the  idea  was  that  a  good  knowledge  of  such  effects  is  necessary  for  a 
correct  interpretation  of  the  results,  in  order  to  distinguish  between  the  real  physical  results  and 
the  artificial  numerical  effects.  Once  we  were  able  to  make  this  distinction,  we  used  our  code  to 
simulate  some  well-known  configurations,  so  as  to  get  more  insight  into  the  physics  of  a  plasma 
sheath. 

In  this  report,  we  first  present  our  physical  and  numerical  models;  we  then  devote  a  sub¬ 
stantial  part  to  the  numerical  methods  used;  we  conclude  with  the  study  of  numerical  and  phy¬ 
sical  effects. 
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2.  MODELS 


2.1.  Physical  Model 

Figure  1  shows  the  configuration  we  will  study  here;  it  can  be  described  as  follows;  the 
half-space  x<0  (we  treat  the  unidimensional  case)  consists  of  a  plasma  reservoir  filled  with 
particles  of  one  or  several  species;  we  assume  that  the  velocity  distributions  of  each  species  in 
the  reservoir  are  Maxwellian  (with  a  possible  drift)  in  the  neighborhood  of  x-0.  A  target  is 
placed  at  x— L>0,  which  absorbs  perfectly  any  incident  particle  (no  secondary  emission). 
Between  x— 0  and  x**£,  we  have  a  vacuum  background,  plus  (possibly)  particles  emitted  from 
the  plasma  reservoir. 

The  virtual  wall  at  x—0  is  assumed  to  be  at  zero  potential;  an  external  potential  difference 
<t>0  can  be  applied  between  x— 0  and  x— I.  <t>0  can  take  any  value  and  can  be  time-dependent. 
Furthermore,  there  can  be  a  net  current. 

We  do  not  add  possible  collisions  between  particles.  No  gravitational  field  is  present. 
Unless  otherwise  specified,  particles  are  electrons.  No  magnetic  field  is  present. 

As  a  consequence  of  the  above  configuration  and  assumptions,  we  can  expect  the  follow* 
ing  general  behavior;  particles  are  moving  about  due  to  their  own  and  applied  fields;  depending 
on  the  charge  density  within  the  vacuum  region,  depending  on  the  externally  applied  potential, 
and  depending  on  its  initial  velocity,  a  particle  contributes  to  the  cathode  (x— 0)  or  anode 
(x—L)  charge  flux,  or  it  remains  between  the  boundaries.  It  is  the  collective  behavior  of  such 
particles,  the  resulting  fields,  currents,  energies,  etc.  that  the  PLOUF  code  allows  us  to  study. 


2.2.  Numerical  Model 

This  paragraph  should  give  the  reader  a  global  idea  of  the  general  structure  of  the  PLOUF 
code,  without  entering  too  many  details  at  this  time.  We  should  point  out  here  that  our  code 
owes  its  structure  and  much  of  its  actual  programming  to  the  ESI  code  written  by  A.  B.  Lang- 
don.  In  this  report,  we  will  assume  that  the  reader  has  some  previous  knowledge  of  this  code  or 
has  access  to  the  relevant  literature. 

2.2.1.  The  Main  Frame 

Since  we  are  studying  time-varying  phenomena  using  a  digital  computer  (as  opposed  to  an 
analog  computer),  some  kind  of  discretization  of  time  has  to  be  introduced:  the  basic  assump¬ 
tion  is  that,  by  knowing  what  happens  at  a  finite  number  of  instants,  we  will  know  the  continu¬ 
ous  evolution  of  the  variables.  Our  code  (as  well  as  ESI)  proceeds  as  follows:  knowing  the  posi¬ 
tions  and  charges  of  the  panicles  at  some  time  t  (and  the  velocities  at  r-<W 2),  we  can  evalu¬ 
ate  the  fields  at  that  time  t,  deduce  the  accelerations,  find  the  new  velocities  at  time  r+dr/2, 
calculate  the  new  positions  at  /+Af,  and  so  on. 

Consequently,  the  main  frame  of  the  code  can  be  very  roughly  described  like  this: 

1  Initialize  (includes:  f— 0;  timestep—br,  duration  of  simulation  terHl;  definition  of  spatial 
grid;  put  panicles  in  their  places  in  phase-space). 

2  Compute  fields  (from  charge  density  p). 

3  Compute  velocities  of  particles. 

4  Compute  positions  of  panicles. 

5  t-f+Ar. 

6  If  t< I,,*,  then  go  to  line  2. 

7  End. 
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2.2.2.  Non-Neutral  Conditions 

The  ESI  code  is  designed  for  the  simulation  of  periodic,  infinite  configurations,  with  no 
net  charge  or  field:  <p>— 0  and  <£>—0.  It  does  not  include  any  source  or  sink  of  particles. 
However,  we  see  from  the  physical  model  described  in  Section  2.1  that  in  our  problem  particles 
can  enter  or  leave  the  sheath  (i.e.,  they  appear  to  be  created  or  destroyed).  In  order  to  allow  for 
a  non-neutral  plasma,  we  let  the  simulation  region  extend  from  x— 0  to  x-2 L,  and  consider 
that  our  true  sheath  is  from  x— 0  to  x— L ,  and  assume  that  the  remaining  half  is  a  fake  sheath. 
In  the  first  region,  we  compute  the  charge  density  p  from  the  particle  positions.  Then,  we  add 
an  inverted  charge  density  in  the  fake  region,  so  that  p (L -f  x)— p (L—x).  This  modification 
means  that  we  have  p(0)— p(2£),  so  that  the  periodicity  requirement  £(x+2I)-£(x)  for  the 
Fourier  transforms  in  ESI  is  satisfied.  Note  that  this  modification  is  equivalent  to  considering 
that,  for  each  charge  at  x—xo  in  the  true  sheath,  there  is  an  image  charge  of  opposite  sign  at 
x—2£— xo  in  the  fake  region. 

Once  we  know  the  charge  density  p(x)  in  the  whole  sheath,  we  can  evaluate  the  fields  for 
0<x<2£,  using  Poisson’s  equation;  then  we  add  the  externally  applied  potential  4>0  between 
x— 0  and  x—I,  as  a  solution  to  the  homogeneous  (Laplace)  equation. 

The  existence  of  a  plasma  reservoir  for  x<0  is  accounted  for  by  the  creation  of  particles 
with  positive  velocities  at  x— 0,  and  by  the  absorption  of  particles  with  negative  velocities  at 
x— 0  and  with  positive  velocities  at  x—I .  (In  fact,  this  absorption  does  not  occur  exactly  at  x— 0 
and  x— I;  see  Section  3.2.2.) 

Note  that  there  are  no  particles  in  the  fake  sheath,  so  that  the  existence  of  this  region 
increases  the  computation  costs  only  in  the  doubling  of  the  region  over  which  the  fields  are 
solved  (in  other  words,  there  is  no  doubling  of  the  number  of  particles). 

With  the  above  features,  PLOUF  is  basically  ready  to  run. 
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3.  THE  PLOUF  CODE 


3.1.  General  Structure 

PLOUF  performs  some  tasks  added  to  ESI,  such  as:  inject  new  particles  in  the  diode, 
compute  currents,  deal  with  the  variable  number  of  particles  in  the  system. 

Consequently,  the  main  frame  of  our  code  looks  like  this: 

1  Initialize  (includes:  /— 0;  timestep—&c,  duration  of  simulation  tend\  definition  of  spatial 
grid;  load  particles  in  phase-space  at  /— 0;  load  entering  panicles  in  their  own  phase- 
space). 

2  Compute  fields  (from  charge  density  p),  and  evaluate  the  displacement  currents  near  the 
boundaries. 

3  Compute  velocities  of  particles. 

4  Compute  positions  of  panicles. 

5  Compute  currents  everywhere  (at  some  selected  time  steps)  and  at  the  boundaries  (at 
each  time  step). 

6  Repack  arrays  representing  particles  (if  applicable) . 

7  Load  entering  particles. 

8  /-t+Ar. 

9  If  then  go  to  line  2. 

10  End. 


3.2.  Initialization 

The  INIT  subroutine  performs  two  tasks.  First,  it  sets  up  the  state  of  the  system  at  the 
stan  of  the  simulation  by  loading  panicles  in  phase-space  according  to  the  relevant  input  param¬ 
eters.  This  set  of  operations  is  performed  exactly  as  in  ESI;  it  should  be  noted  that  the  initiali¬ 
zation  stage  is  not  of  primary  concern  to  us,  since  our  problem  is  essentially  a  boundary-value 
problem  (not  an  initial- value  problem  like  ESI). 

The  second  task  is  to  load  the  new,  entering  panicles  in  phase-space.  Specifying  how 
many  panicles  will  enter  the  system  during  each  time  step,  we  set  up  (in  SETVNU)  a  discrete 
distribution  of  velocities  for  these  panicles.  Then,  we  compute  (in  SETXNU)  the  corresponding 
initial  positions,  which  cannot  be  merely  x— 0,  due  to  the  method  we  use  for  computing  the 
currents  at  the  boundaries. 


3.2.1.  SETVNU 

According  to  our  physical  model  (see  section  2.1),  particles  in  the  reservoir  x<0  have  a 
(possibly  drifting)  Maxwellian  distribution  of  velocities.  Under  this  assumption,  the  purpose  of 
SETVNU  is  to  compute,  for  each  species  independently,  the  velocities  of  the  particles  entering 
the  system  at  each  time  step,  given  the  following  parameters: 

nenter  the  number  of  particles  n^,  in  the  distribution; 
vinth  the  thermal  velocity  v,*  of  the  distribution; 
vinbm  the  drift  velocity  vb  of  the  distribution; 
vinmin  the  smallest  velocity  vmn  in  the  distribution; 

vinmax  the  largest  velocity  in  the  distribution.  (These  last  two  parameters  allow  retaining 
only  a  "slice"  of  the  whole  Maxwellian.) 

Let  us  see  how  these  nenter  velocities  are  computed.  In  the  cold  case,  i.e.  vtH-0,  the 
result  is  trivial: 
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In  the  warm  case,  the  task  is  more  complex.  Lets  consider  that  we  have  the  following  normal¬ 
ized  Maxwellian  distribution  in  the  reservoir  x<0: 

-V2 

J — e  dv„w 

virv,A 

Then  the  distribution  function  for  emission  (i.e.  the  probability  to  observe  during  ±i  a  particle 
passing  through  jc-0  with  a  velocity  between  vMW  and  v„„+t with  v>0)  is: 

g(v^)dv^,-^l"—varwe  ^  dv„w  (1) 

where  A  is  a  normalization  constant  defined  by: 


ft  new 


Changing  variables,  x— we  obtain: 


Now  that  g(v)dv  is  specified,  SETVNU  executes  the  following  algorithm,  which  constructs  a 
set  of  nenter  velocities  satisfying  Equation  1: 


1 

2 


3 


Initialize:  /-I 

Compute  the  integral  of  the  distribution  function  g(v)dv  from  zero  to  the  velocity  of  the 
slowest  particle  to  be  injected;  this  will  give  the  value  of  the  integral  corresponding  to  the 
smallest  velocity  in  the  distribution  function  for  entering  particles: 


f  min” 


I 

Rnew  1 


Compute  the  value  of  the  integral  of  g{v)dv  which  corresponds  to  the  next  higher  velo 
city  in  the  distribution: 


2/ 


4  Using  Newton’s  method,  compute  the  velocity  corresponding  to  the  value  of  the  integral 

v 

that  we  just  found  in  Step  3;  that  is,  find  v  such  that  Jf(»w)dvw-/ 

*o 

5  Prepare  next  step:  /“/+ 1 

6  If  /  then  go  to  line  3 

7  End 

In  Step  4  above,  using  Newton's  method  implies  having  some  knowledge  of  an  initial 
approximation  v0 ( / )  to  the  solution  >*,(/).  The  safest  way  to  go,  when  /-l,  is  to  start  with 

Vo 0)“ - 2 -  ’ 

i.e.  the  velocity  for  which  g(vww)  is  maximum.  This  way,  Newton’s  method  converges  even  for 
sharply  peaked  distributions  (without  becoming  too  costly  in  terms  of  computing  time).  Subse¬ 
quently,  for  />!,  we  just  take  v0(/)«  vw(/-l). 
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3.2.2.  SETXNU 

The  SETXNU  subroutine  computes  the  positions  of  the  entering  panicles;  these  positions 
are  computed  according  to  the  following  three  criteria: 

*  all  new  panicles  must  pass  x— 0  regularly  in  time;  this  means  for  instance  that  the  fth  panicle 
of  one  group  must  enter  the  diode  at  time  t+Ati(i— D/nenter]  (assuming  that  we  have  one  sin¬ 
gle  species).  Note:  we  call  "group*  a  set  ot  nenter  panicles  of  one  or  several  species,  due  to 
enter  the  diode  during  one  time  step  of  the  simulation. 

*  all  species  must  be  treated  simultaneously  within  one  group;  e.g.:  if  we  have  two  species  with 
nenter Inenter  ■*,  then  we  will  first  load  two  panicles  of  species  1,  then  one  panicle  of  species 

2,  then  two  panicles  of  species  l . until  nenter  |+ nenter 2  particles  have  been  loaded  for  the 

given  time  step. 

*  no  panicle  may  be  loaded  at  x>— Ax/2;  this  is  due  to  our  method  of  current  measurements. 

On  the  left-hand  side  of  Figure  2,  we  show  the  positions  of  the  panicles  when  they  are 
loaded  in  phase-space,  assuming  (for  the  sake  of  simplicity)  that  the  distribution  function  for 
emission  is  such  that  the  free-streaming  panicles  form  a  straight  line  in  phase-space.  Under  this 
assumption,  the  distribution  function  tells  us  that  the  panicles  are  to  be  loaded  along  a  straight 
line  in  phase-space.  However,  it  does  not  tell  us  either  the  slope  or  the  position  of  this  line, 
which  are  determined  as  follows: 

*  the  position  of  the  line  (its  offset  in  the  x-direction  once  the  slope  is  known)  must  be  such 
that  the  slowest  panicle  is  loaded  at  a  position  xOo)**— Ax/ 2,  where  r0  is  the  loading  time; 

*  calling  n  the  number  of  time  steps  needed  for  the  slowest  panicle  to  enter  the  sheath  from  its 
loading  position,  then  the  slope  of  the  straight  line  must  ensure  that  exactly  one  panicle  passes 
xH)  at  time  t0+nAt+At/nenter,  exactly  one  panicle  passes  x— 0  at  time 
t0+nAt+2At/nenter,  ....  and  exactly  one  panicle  passes  x-0  at  time 
fo+  n  A  t 4- nenter  A  //  nenter . 

The  curve  on  the  right-hand  side  of  Ftgure  2  (the  curve  located  inside  the  diode  region) 
shows  the  positions  of  the  particles  in  phase-space  at  time  r0+(n+l)Af,  i.e.  at  the  end  of  the 
time  step  during  which  they  enter  the  sheath. 

Figure  3  is  similar  to  Figure  2,  except  that  it  corresponds  to  a  Maxwellian  distribution 
function  in  the  reservoir. 


3.3.  HELDS 

The  subroutine  FIELDS  has  been  modified  so  as  to  construct  the  inverted  charge  density 
in  the  fake  sheath,  and  to  add  the  externally  applied  potential.  For  further  explanations,  please 
read  Section  2.2.2. 

This  subroutine  also  computes  the  displacement  currents  near  the  boundaries,  i.e. 
(Ax)  and  Jd,v(L-Ax): 

,  /...x  E(Ax,t+At)-E(Ax,t) 

Jdisp  l  Ax )  - - 

,  ,,  .x  E(L-Ax,t+At)—E(L-Ax,t) 

/«*  (I -Ax)- - — - 

Please  read  Section  3.5.2  to  see  why  we  compute  /,,*( Ax)  and  Jdlsp(L-Ax)  instead  of  Jdnp (0) 
and  /**(/). 
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3.4.  ACCEL 


The  ACCEL  subroutine  computes  the  new  velocities  at  time  r+A// 2. 

Sections  3.5. 1.2  and  3. 5. 1.3  explain  how  the  currents  are  measured.  A  consequence  of 
our  method  of  loading  new  particles  is  that  they  are  not  exactly  created  and  destroyed  at  x— 0 
and  x They  are  in  fact  created  at  some  x<-Ax/2  and  destroyed  at  x— lx/2 
andx-Z.  +Ax/2.  Of  course,  this  slight  modification  with  respect  to  our  numerical  model  must 
not  have  any  consequence  on  the  behavior  of  physical  quantities  between  x— 0  and  x— L.  The 
ACCEL  subroutine  takes  care  of  this  by  keeping  constant  the  velocities  of  the  particles  which 
are  at  positions  xo  such  that  -Ax/2^xo<0  or  L  <xo<L+Ax/2. 

In  other  words,  particles  free  stream  when  they  are  n. .  strictly  within  the  diode  (i.e.  when 
their  positions  x  are  not  such  that  0<x<L).  Note  that  this  is  an  approximation,  since  we  have 
finite-size  particles,  which  should  be  accelerated  proportionally  to  that  fraction  of  the  particle 
which  is  strictly  within  the  diode. 


Furthermore,  ACCEL  computes  the  drift  and  kinetic  energies  so  that  the  contribution  of 
a  particle  is  multiplied  by  a  factor  which  gives  the  percentage  of  the  finite-size  particle  residing 
inside  the  diode: 


part(x) 


I+Ax 

L 

,  1 

max 

2 

2  ”*| 

,  o 

This  accounts  for  the  gradual  entrance  of  a  finite-size  particle  in  the  diode,  and  it  decreases  the 
noise  (in  the  energy  calculations)  due  to  particles  exiting  and  entering  the  diode. 


3.5.  MOVE 

The  basic  way  of  computing  the  new  positions  of  the  particles  at  each  time  step  is  trivial: 
x(r+A/)—x(/)+yr(r+Ar/2)Ar.  The  mover  computes  the  new  positions  of  all  particles  accord¬ 
ing  to  this  relation,  as  in  ESI. 

In  our  code,  MOVE  has  also  a  second  function,  which  is  to  measure  the  charge  fluxes 
(also  called  ballistic  currents). 


3.5.1.  Computation  of  the  Charge  Fluxes 

There  are  different  means  of  obtaining  the  ballistic  currents.  One  method  would  be  to 
compute  j—pv  knowing  v  at  time  r+Ar/2  and  evaluating  p  at  the  same  time,  using  pit)  and 
p(r+Ar).  This  method  requires  that  we  define  some  interface  (some  interpolation  scheme) 
between  the  velocity,  which  is  a  particle  quantity,  and  the  charge  density,  which  is  a  grid  quan¬ 
tity;  furthermore,  we  will  see  in  Section  3. 5. 1.1  that  this  current  does  not  satisfy  KirchholTs  law 
within  any  volume  bounded  by  any  pair  of  grid  planes. 

The  method  implemented  in  PLOUF  measures  i—dq/dt,  i.e.  the  charge  which  has  gone 
through  the  grid  plane  (s)  on  which  we  want  to  measure  the  charge  flux. 

At  this  point,  please  allow  us  to  briefly  recall  some  basic  notions  about  particle  shaping. 

In  order  to  evaluate  the  particle  (and  charge)  density  on  the  discrete  spatial  grid  from  the 
continuous  particle  positions,  we  infer  that  the  particles  have  some  finite  size.  We  call  particle 
shaping  (or  particle  weighting  )  the  process  in  which  we  define  (on  the  grid)  a  particle  density 
n(x)  for  each  particle.  Typically,  the  extent  of  a  particle  (i.e.  the  length  of  the  interval  on 
which  n  (x)  s*0)  is  one  or  two  grid  spaces. 
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3.5.1. 1.  The  >—pv  Method 

Now,  let  us  examine  the  relation  between  the  particle  shaping  and  the  measurement  of 
the  ballistic  current  via  the  first  method  O—pv). 

We  need  to  evaluate  p,+i,/2  on  the  basis  of  p,  and  p,+i/.  We  can  take  p,*i(/ 2  to  be  the 
time-average  of  the  shapes  p,  and  p,+i(  (see  Figures  4  and  5): 


Pr+\t/2^Xj  t  Xt,  ^(•*■<1/)  —  P 


x,+x,+x, 

2 


We  can  also  define  p  ,4*1/2  as  the  space-average  of  p,  and  p,+i(  (see  Figures  6  and  7): 

,  _  p(X,,xr)  +  p(Xj,x.^Xl) 

The  advantage  in  using  either  one  of  these  definitions  is  that  we  will  measure  instantane- 
ous  charge  fluxes  with  a  high  accuracy.  In  fact,  the  time-averaging  for  p(xr4.d{/2)  would  be 
required  to  insure  self-consistency  if  we  were  to  re-use  these  ballistic  currents  at  some  other 
place  in  the  simulation  (e.g.  to  compute  new  fields);  in  PLOUF,  the  currents  are  needed  only 
for  diagnostics. 

The  bad  point  is  that  this  method  does  not  provide  for  charge  conservation  within  every 
volume  bounded  by  an  arbitrary  choice  of  two  planes.  Indeed,  looking  at  what  happens  for  a 
high-velocity  particle  (Figures  2a  and  3a),  we  see  that  there  may  be  grid  planes  with  no  charge 
fiux,  even  though  in  reality  some  charge  has  gone  through  it  in  one  time  step.  At.  In  other 
words,  KirchhofFs  law  is  not  satisfied  with  this  model,  and  there  is  no  guarantee  (especially  for 
randvx—Q)  that  this  adverse  effect  will  become  statistically  negligible.  Furthermore,  some  of 
the  dynamics  would  probably  be  lost  in  the  interface  between  charge  densities  and  particle  velo¬ 
cities. 


3.5.I.2.  The  i-dq/dt  Method 

The  second  method  settles  the  problem  of  charge  conservation,  at  the  expense  of  some 
loss  of  information  about  current  dynamics.  This  method  computes  the  total  charge  flow 
through  each  grid  plane  during  one  time  step: 

.  __  Aq  lf;.# 

'mtasurrd  —  J  l,ral 

where  is  defined  on  grid  planes  only: 

*measurtd  —  1 measurrd  0  Ax) 

and  is  the  instantaneous  charge  flux  at  any  given  instant  between  t  and  t+At. 

In  order  to  achieve  this  evaluation  of  the  current,  we  must  use  a  current  shape  as  defined 
by  Figures  8  and  9.  Using  this  shape  and  taking  its  values  on  the  grid  points,  we  have  directly 
the  values  of  the  currents  at  these  points. 


3.5.I.3.  Current  Measurements 

After  having  moved  a  particle,  MOVE  computes  the  ballistic  currents  associated  with  this 
particle;  it  then  gives  two  results: 


a  values  i{Xj,t)  on  all  grid  planes,  at  some  selected  time  steps; 

•  boundary  charge  fluxes:  itmmrU), 

Like  the  accumulation  of  charge  density  ,  the  computation  of  the  current  shapes  is  basi¬ 
cally  non-vectorizable,  thus  making  MOVE  quite  expensive  to  use.  However,  we  do  not  really 


Figure  4  (Section  3.5. 1.1):  Time-avenging  of  two  particle  shapes  (fast  particle). 
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need  to  know  the  currents  everywhere  at  each  time  step.  (In  fact,  the  resulting  plot  would  be  so 
overfull  of  information  that  it  would  be  almost  entirely  black!)  So,  these  computations  are  per¬ 
formed  only  once  every  nth/S  time  steps,  where  nth  is  the  number  of  time  steps  between  his¬ 
tory  plots;  in  other  words,  each  three-dimensional  plot  of  the  ballistic  currents  will  have  six 
curves  i(x). 

On  the  other  hand,  we  also  want  history  plots  of  the  boundary  charge  fluxes  (emitter, 
cathode  and  anode  fluxes).  Here  it  would  be  insufficient  just  to  compute  the  ballistic  currents 
every  (nth/5)th  time  step,  but  since  these  currents  are  related  to  particles  passing  isolated  grid 
planes,  the  problem  of  constructing  the  current  shapes  is  no  more  relevant,  and  the  computa¬ 
tion  is  now  vectorizable.  Consequently,  the  boundary  charge  fluxes  are  evaluated  independently 
(from  the  ballistic  currents  in  the  whole  diode)  in  a  second  part  of  the  subroutine. 


3.S.2.  Computation  of  the  Displacement  Currents 

Displacement  currents  are  in  fact  computed  in  FIELDS,  but  we  first  needed  to  introduce 
the  current  shapes  (Section  3.5. 1.2),  in  order  to  explain  why  these  variations  of  the  electric 
field  are  evaluated  at  x-Ax  and  x— L-Ax  instead  of  x-0  and  x-L. 

Let  us  look  at  an  example.  We  consider  a  single  particle  just  about  to  enter  the  diode  at 
x— 0;  we  assume  that  Lx/2,  and  that  — Ajr/2<jc,0+ix,<0.  Then,  the  electric  field  due  to 
this  particle  at  times  t0  and  r0+Ar  is  zero  at  both  times  (Figure  10),  so  that 

,  £(/o+Af)— £(fo) 

d,sp  At 

However,  as  one  can  see  in  Figure  11,  the  charge  flux  that  we  obtain  by  using  the  method 
described  in  Section  3.5.2. 1  is  nonzero  at  the  emission  plane. 

So,  in  order  to  reduce  this  discrepancy  between  the  ballistic  and  displacement  currents,  all 
’boundary”  currents  are  computed  one  grid  cell  inside  the  diode  (not  on  the  boundary  them¬ 
selves).  Consequently,  we  call: 

-  emitter  charge  flux:  i—dq/dt  due  to  particles  passing  X^—Ax  with  vx>0; 

-  anode  charge  flux:  i—dq/dt  due  to  particles  passing  X hG—L-Ax  with  vr>0; 

-  cathode  charge  flux:  i—dq/dt  due  to  particles  passing  X:-Ax  with  yT<0; 

-  currents  at  emitter  ballistic,  displacement  and  total  currents  at  X:-Ax\ 

•  currents  at  target:  ballistic,  displacement  and  total  currents  at  XvC—L- Ax. 


3.6.  REPACK 

The  number  of  particles  in  the  system  is  not  constant  in  time.  However,  some  arrays  in 
the  code  have  sizes  proportional  to  the  total  number  of  simulated  particles.  Hence  we  must 
define  the  sizes  of  these  arrays  in  such  a  way  that  we  never  attempt  to  toad  particles  beyond  the 
limits  of  the  arrays.  On  one  hand,  a  straightforward  solution  to  this  question  is  to  repack  the 
arrays  from  time  to  time,  and  to  expand  them  whenever  we  try  to  go  above  their  upper  bounds. 
On  the  other  hand,  we  could  declare  their  sizes  to  be  much  targer  than  the  initial  number  of 
particles,  adding  new  particles  at  the  end  of  the  area  (of  the  array)  already  occupied  by  other 
particles.  The  first  solution  is  very  costly,  and  the  second  one  sets  a  very  stringent  limit  on  the 
number  of  time  steps  allowed  in  one  run. 

The  REPACK  subroutine  executes  the  following  compromise:  the  size  of  the  arrays  are 
chosen  equal  to  (or  greater  than)  the  maximum  number  of  particles  that  will  ever  be  simultane¬ 
ously  in  the  system.  Then,  the  contents  of  the  arrays  are  repacked  every  ipackth  time  step, 
where  ipack  is  given  as  an  input  parameter  (a  typical  value  is  ipack~ 8). 
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3.7.  ENTER 

The  role  of  the  ENTER  subroutine  is  to  load  new  particles  into  the  system  during  the 
simulation.  The  velocities  of  the  entering  particles  are  those  computed  in  SETVNU,  and  their 
initial  positions  have  been  evaluated  in  SETXNU. 

With  this  method,  we  see  that  entering  particles  have  the  same  nenter  initial  velocities  at 
each  time  step.  This  does  not  take  into  account  the  noise  effects  due  to  the  existence  of  a 
semi-infinite  plasma  reservoir  in  the  half-space  x<0,  and  it  results  in  a  multibeaming  injection, 
rather  than  a  smooth  one  (averaged  over  time).  This  problem  has  been  reduced  by  the  intro¬ 
duction  of  a  small  noise  (different  at  each  time  step)  in  the  distribution.  This  noise  is  intro¬ 
duced  in  the  form  of  a  randomization  added  to  the  velocities  given  by  SETVNU;  the  amplitude 
of  this  noise  (i.e.  the  amount  of  randomization  added)  is  given  by  the  input  parameter  randvx 
(see  Section  4.1.1). 


3.8.  Input  Parameters 

PLOUF  allows  the  user  to  control  a  large  number  of  physical  and  simulation-related  vari¬ 
ables.  We  describe  here  briefly  the  most  important  of  these  parameters. 

Definition  of  the  system: 
dt  length  of  the  time  step; 

/  physical  length  of  the  true  sheath; 

ng  number  of  grid  cells; 

nt  number  of  time  steps  during  the  whole  simulation; 
phiO  externally  applied  potential,  which  can  be  time-dependent  (see  Figure  12); 
phiQ(lj)  is  the  y'th  given  value  of  <t>0; 
phiQ(2j )  is  the  time  f|  when  <t>0sph/0(l  J)  must  be  reached; 
phiOOj)  defines  the  evolution  of  <M*)  between  times  t0"mphiO(2j— 1)  and 

fi-ph/0(2j). 

randvx  percentage  of  noise  added  to  to  the  ordered  velocities. 

Definition  of  each  species: 

nenter  number  of  new  particles  at  each  time  step; 

qm  charge-over-mass  ratio; 

vmbim  drift  velocity  of  the  Maxwellian  (in  the  reservoir); 
vinmax  maximum  velocity  of  entering  particles; 
vinmin  minimum  velocity  of  entering  particles; 

W nth  thermal  velocity  of  entering  particles; 

w p  plasma  frequency 


3.9.  Diagnostics 

PLOUF  gives  essentially  two  kinds  of  graphic  outputs:  snapshots,  which  give  the  states  of 
the  sheath  at  some  given  instants,  and  history  plots,  which  describe  the  time-evolution  of 
different  variables. 


Figure  13  (Section  4.1.1):  Graphical  explanation  of  the  definition  of  randvx. 
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3.9.1.  Snapshots 

The  snapshots  in  PLOUF  are  the  same  as  in  ESI.  They  are  mainly: 

-  charge  density  p  as  a  function  of  position  x; 

•  electric  potential  $(x); 

•  electric  field  £(x); 

•  positions  of  particles  in  phase-space  (x,vx)\ 

•  velocity  distribution  /(vt). 


I 
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3.9.2.  History  Plots 

We  list  below  the  diagnostics  given  by  the  history  plots  in  our  code  (the  plots  marked 
with  an  “**  are  new  features  not  to  be  found  in  ESI): 


*  externally  applied  potential  <t>0({):  this  displays  the  time-evolution  of  4>0(t)  as  prescribed  by 
the  input  array  phiQ. 

*  plasma  frequency  a»„(r):  gives  the  plasma  frequency  averaged  over  the  whole  diode,  as  a 
function  of  time; 

*  potential  minimum  $min0 r,/):  shows  the  time-evolution  of  the  position  x  and  the  value  4>min 
of  the  minimum  of  the  potential.  (This  minimum  is  interpolated  between  three  neighboring 
grid  points.) 

*  charge  fluxes  i(xJtt):  shows  the  time-evolution  of  i(xj). 

*  emitter  charge  flux:  plots  the  ballistic  current  due  to  particles  passing  through  the  grid  plane  at 
Xi^^x  with  positive  velocities. 

*  anode  charge  flux:  plots  the  ballistic  current  due  to  particle  passing  through  the  grid  plane  at 
X\q  L  —  Ajc 

*  cathode  charge  flux:  plots  the  ballistic  current  due  to  particles  passing  through  the  grid  plane 
at  Xf“\x  with  negative  velocities. 

*  charge  fluxes  at  boundaries:  groups  together  the  previous  three  plots,  on  the  same  scale. 
Note:  all  the  charge  flux  plots,  except  i(Xj,t),  can  be  displayed  separately  for  each  particle 
species  and/or  for  all  species  together. 

*  currents  at  emitter:  plots  the  ballistic,  displacement  and  total  currents  at  T2-4jc. 

*  currents  at  target:  plots  the  ballistic,  displacement  and  total  currents  at  XNG—L- Ax. 


-  field  energy:  £> (/)  — 

*  2 

. .  „  .  .  _  m,  v‘ 

-  kinetic  energy:  £*(r)  -  £ — - — . 

'  .7  2 

-  drift  energy:  £0 (t)  -  £  -  ' 

/  **  _ 

ff\  (  y  ^— *) 

-  thermal  energy:  £>(;)  -  £ — : — — — 

1  ~ 

-  total  energy:  £j(/)  -  £>  +  £*. 
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4.  RESULTS 


4.1.  Numerical  Results 

4.1.1.  Initial  Velocity  Distribution  and  Noise 

As  mentioned  earlier  (in  section  3.2.1),  particles  entering  the  sheath  at  x—0  have  their 
velocities  chosen  so  as  to  simulate  a  (possibly  drifting)  Maxwellian  distribution  for  x<0,  with 
possible  randomization.  The  amount  of  added  noise  is  controlled  by  the  input  parameter 
randvx,  the  effects  of  which  we  want  to  discuss  now. 

randvx  is  what  we  define  as  the  percentage  of  noise  added  to  the  regular  Maxwellian.  If 
vx^. .Maekssd—l)  and  wc 0)  are  two  neighboring  velocities  of  the  noiseless  Maxwellian 
defined  in  part  3.2.1,  then  randvx  allows  choice  of  vxmsy(i)  (the  initial  velocity  actually 
assigned  to  the  particle)  randomly  between  these  values  using  the  algorithm 

VXfm.mmy  (()  m  vxntw .noiitka  0)"”[  VX»ew,,»urlfu  0)~  VX^,_ imstltu  U“l)]randVX  . 

Figure  13  shows  that  randvx  defines,  for  each  aOT,«wiw('),  the  range  of  velocities  within 
which  vx _ „.(/)  will  be  randomly  chosen. 

Of  course,  such  a  noise  will  increase  the  noise  levels  of  other  quantities  measured  within 
the  diode.  However,  if  randvx  is  kept  reasonably  low  (randvx*  1  is  OK  in  many  cases),  the  phy¬ 
sical  results  are  generally  unchanged  qualitatively.  The  disadvantage  of  adding  noise  is  obvious 
(adding  noise  is  noisy!);  it  may  however  be  desirable  to  do  so  in  order  to  suppress  nonphysical 
diagnostics. 

Figure  14  shows  that  too  small  a  value  of  randvx  yields  a  multistreaming  phenomenon  in 
phase-space:  each  line  corresponds  to  one  velocity  in  the  noiseless  Maxwellian. 

Figure  15  shows  that  the  addition  of  a  small  noise  yields  a  much  more  physically  plausible 
diagnostic:  particles  are  indeed  confined  to  some  regions  in  phase-space,  but  vx(x)  is  an 
infinitely-multiple-valued  function  of  x. 

Figure  16  shows  how  rapidly  (and  beautifully)  the  equilibrium  current  is  reached  (time 
r— 0  is  in  the  back,  r— 5  is  in  the  front). 

Figure  17  shows  that  the  equilibrium  is  still  reached  very  rapidly  (see  section  4.2.2),  but 
that  it  remains  very  noisy  when  randvx  ^0. 

The  above  results  show  that  it  may  be  desirable  to  add  noise  for  some  purposes,  whereas 
it  just  gives  dirtier  results  in  other  cases.  We  have  also  found  multistreaming  in  phase-space  to 
be  a  nonphysical  effect. 


4.1.2.  Oscillations  of  the  Potential  Minimum  (Driftless  Case) 

Previous  results  showed  that,  for  a  Maxwellian  with  no  beam  component  (i.e.  nondrifting; 
in  the  reservoir,  the  equilibrium  solution  is  non-oscillatory. 

However  our  own  results  almost  always  display  some  kind  of  time-dependence  of  the 
equilibrium  state,  which  can  be  seen  in  the  variation  (in  amplitude  and  position)  of  the  poten¬ 
tial  minimum  (Figure  18).  Even  though  this  variation  is  quite  periodic  in  some  instances, 
a  closer  look  at  it  shows  that  it  is  not  a  physical  phenomenon,  but  a  numerical  effect  due  to  the 
lack  of  resolution  inherent  to  our  simulations. 


Looking  at  the  motion  of  in  (x,<t>) -space,  we  have  noticed  the  following  effects: 

*  The  oscillation  of  $mtn  occurs  around  the  values  obtained  on  the  basis  of  Langmuir’s  findings. 

*  The  x-variation  of  spans  only  a  few  inter-particle  spaces,  generally  two  or  three.  This 
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means  that,  with  the  particle  densities  we  had  (corresponding  to  simulation  runs  with  10,000  to 
100,000  particles),  the  position  of  <t>min  was  basically  depending  on  information  furnished  by 
only  three  or  four  neighboring  particles.  This  is  a  very  crude  way  of  finding  the  true  potential 
minimum. 

*  The  amplitude  of  the  variation  of  4>min  decreases  when  the  number  of  particles  increases,  i.e. 
when  the  particle  density  around  x#  increases  (Figure  19).  This  is  in  part  a  consequence  of 
the  previous  observation. 

*  The  behavior  of  1>min  is  not  very  sensitive  to  the  number  of  grid  cells. 

The  above  results  show  that  the  variations  of  4>min  and  x*^  are  essentially  a  consequence 
of  the  coarseness  of  the  simulation.  Increasing  the  number  of  particles  and  decreasing  At,  we 
tend  asymptotically  towards  a  time-independent  state.  Thus  these  oscillations  of  the  potential 
minimum  are  a  numerical  effect. 

Note:  this  oscillation  has  never  been  found  to  give  rise  to  (numerical)  particle  trapping. 
Indeed,  a  necessary  condition  for  a  particle  to  be  trapped  due  to  the  motion  of  <t>mm  is  that  the 
trajectory  of  <bmM(x,r)  has  two  maxima  in  x*  -space  (Figure  20);  but  in  our  case,  this 

(TUB 

trajectory  has  only  one  maximum  (Figure  21),  thus  forbidding  any  trapping. 


4.1.3.  Time-Dependence  of  the  Currents 
4.1.3. 1.  Current  Noise 

We  have  seen  in  Section  4.1.1  that  the  currents  become  noisier  when  randvx  is  increased. 
Our  datas  show  that  the  emitter,  anode  and  cathode  current  noises  have  an  essentially-linear 
dependence  on  randvx. 

We  have  also  compared  the  noises  of  these  three  boundary  charge  fluxes  measured  with 
the  i—dq/dt  method  described  in  part  3.4.2.2,  with  their  values  measured  with  a  zeroth-order 
method  (i.e.  one  in  which  the  ballistic  currents  are  computed  assuming  point-particles).  We 
found  our  own  method  to  decrease  the  noises  by  factors  of  2  to  100;  the  higher  the  particle 
velocities  at  the  boundaries,  the  lower  our  improvement  factor. 


4.I.3.2.  Current  Oscillatioos 

One  can  easily  see  from  Sections  3.5. 1.2  and  3.5.2  that  the  measured  currents  are  sensi¬ 
tive  to  the  grid  spacing  Ax  and  the  length  of  the  time  step  At,  as  well  as  to  the  particle  density 
and  the  particle  velocities  at  the  plane  where  we  measure  the  currents. 

So,  whenever  some  kind  of  oscillatory  solution  is  found  from  the  simulations,  chances  are 
that  it  is  a  numerical  phenomenon,  not  a  physical  one,  if  any  one  of  the  following  conditions  is 
met: 

•  if  A t  is  doubled,  then  the  frequency  of  the  oscillation  is  halved,  and  its  amplitude  is  doubled; 

•  if  Ax  is  doubled,  then  the  frequency  of  the  oscillation  is  doubled,  and  its  amplitude  is  halved; 
>  if  the  particle  density  is  doubled,  then  the  frequency  of  the  oscillation  is  doubled,  and  its 
amplitude  is  halved; 

•  if  the  average  particle  velocity  around  the  plane  used  to  measure  the  currents  is  doubled,  then 
the  frequency  of  the  oscillation  is  doubled,  and  its  amplitude  is  halved. 

If  all  the  above  tests  are  negative,  then  the  oscillation  is  most  likely  some  kind  of  physical 
effect. 

Note:  the  reader  may  have  noticed  that  the  oscillation  of  the  potential  minimum,  as 
described  in  Section  4.1.2,  would  have  failed  the  test  for  Conditions  1  and  3. 


-29- 


i  Figure  19  (Section  4.1.2):  Amplitude  A  of  the  motion  of  the  potential  minimum,  as  a  func- 

(  tion  of  the  particle  density. 
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4.1.4.  Cost  of  the  Pul-Cells 

We  call  pad-cells  these  cells  of  length  Ajc/2  that  we  add  at  both  ends  of  the  true  sheath  in 
order  to  have  correct  current  diagnostics  at  the  boundaries  (see  Section  3.5. 1.2). 

The  cost  of  adding  such  cells  is  reflected  into  an  increase  of  the  totai  number  of  particles. 
Since  Ix—L/ng,  this  cost  increase  will  obviously  depend  on  ng.  We  find  this  increase  to  range 
from  5%  (for  ng— 16384,  with  L— ir)  to  33%  (for  ng—64)  of  the  initial  computation  cost. 

The  charge  flux  measurement  method  described  in  section  3. 5. 1.2  causes  a  dependence  of 
the  current  noise  on  the  grid  spacing:  increasing  the  number  of  ceils  increases  the  noise  level. 
We  found  ng— 512  to  be  the  best  compromise;  in  this  case,  the  pad-cells  represent  an  increase 
of  about  7%  of  the  computation  costs. 


4.2.  Physical  Results 

4.2.1.  Vacuum  Diode 

The  well-known  Child’s  law  has  been  successfully  verified  (cold  injection,  £(0)— 0.): 

4>0c)  -*o(*/I)4/3. 

Furthermore,  we  simulated  the  classical,  warm,  space-charge-limited  diode,  with  half- 
Maxwellian  injection,  and  compared  our  particle  code  diagnostics  with  analytic  results  obtained 
by  William  Lawson,  based  on  steady-state  diode  theory  as  developed  originally  by  Langmuir. 
Lawson  used  a  simple  program  package  (not  a  simulation  code)  to  find  the  potential  and  other 
relevant  quantities  inside  a  unidimensional  (plane)  diode  for  a  limited  set  of  parameters;  for  his 
analysis,  he  had  to  follow  Langmuir  and  assume  that  the  potential  had  a  minimum  somewhere 
inside  the  diode. 

We  computed  the  positions  and  values  of  the  potential  minimum  for  a  set  of 
different  input  values.  Our  results  agreed  quite  accurately  with  Lawson’s  findings.  The  potential 
minimums  obtained  with  our  simulations  were  all  less  than  LI  100  apart  from  the  positions 
predicted  by  Lawson;  note  that  this  was  obtained  with  ng— 64,  i.e.  one  grid  cell  was  more  than 
U 100. 

Figure  22  shows  the  potential  4>(x)  in  a  typical  space-charge-limited  diode,  with 
d>U)  -0. 

4.2.2.  Reaching  Equilibrium 

We  measured  how  rapidly  equilibrium  was  reached  for  time-independent  4>0,  and  tested 
the  uniqueness  (with  respect  to  d>0)  of  such  an  equilibrium. 

Let  us  first  define  our  plasma  frequency: 

Note  that  the  time  scales  on  all  our  (single-species)  plots  is  not  in  unit  of  but  in  units  of 
<Vo: 


A  rule  of  thumb  for  the  plots  featured  in  this  report  is:  =8.0w? 0. 

We  found  that,  in  all  time-independent  cases  that  we  tested,  equilibrium  is  quite  rapidly 
obtained:  using  It— 0.01  (At  is  the  length  of  one  time  step  in  units  of  equilibrium  was 
reached  as  early  as  three  plasma  periods  T„ 3«7.«r The  t‘me  required  to  reach  equilibrium 
increases  with  vmbim/vinth  (see  Section  4.2.3). 
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Figure  22  (Section  4.2.1):  Space-charge-limited  vacuum  diode. 
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Figure  23  shows  the  time-variation  of  the  total  energy  during  the  settlement  of  equili¬ 
brium.  At  r— 2.5  r,.o,  <*>o  is  abruptly  changed  to  show  the  settlement  of  a  second  equilibrium. 

Figure  24  shows  the  motion  of  the  potential  minimum.  At  r— 0,  0  due  to  our 

initialization.  4>mm  then  moves  towards  its  equilibrium  position;  we  have  shown  in  Section  4.1.2 
that  in  fact  4>min  is  numerically  time-dependent  and  circles  in  (x,<t>)-space  around  the  value 
given  by  classical  theory. 

Figure  16  plots  one  curve  i(x)  every  0.5  Tp  #5.  The  curve  in  the  back  represents  the 
current  from  r—0  to  r— dr,  the  curve  immediately  in  front  of  it  shows  the  current  from 
r— 0.5  Tp  a  to  t— 0.5  T^.o+il/.  One  sees  that  all  curves  in  front  of  these  two  (i.e.  obtained  at  later 
times)  already  correspond  to  an  equilibrium  state. 

We  also  tested  the  uniqueness  of  equilibrium  with  respect  to  4>0  (i.e.:  are  there  two  or 
more  possible  equilibria  for  some  values  of  <t>o?).  We  first  reached  a  stable  equilibrium,  then  let 
4>o  vary  slowly  from  its  equilibrium  value  ,  then  went  back  to  ;  we  repeated  this  pro¬ 
cess  for  several  values  of  <t>0  (in  particular  values  such  that  w*.-0  for  4>o-4>lW.  -64*  and 
for  +54>).  We  did  not  find  any  hysteresis-like  phenomenon  (fortunately!);  the 

success  of  this  test  supports  our  claim  that  the  way  the  phase-space  is  filled  at  r—0  is  intrinsi¬ 
cally  irrelevant;  at  worst,  it  could  lengthen  the  time  it  takes  to  reach  equilibrium. 


4.2*3.  Temperature  Effects 

Regardless  of  the  temperature  T,  one  has  always  either  one  of  the  following  two  cases: 

•  all  particles  reach  the  anode:  the  anode  charge  flux  equals  the  emitter  charge  flux  (at  equili¬ 
brium),  the  potential  minimum  is  virtually  motionless,  etc.  The  solution  is  completely  time- 
independent. 

•  some  or  all  of  the  particles  contribute  to  the  cathode  charge  flux:  this  is  the  case  that  we  want 
to  study  below. 

When  'r„,Ao*?*0,  our  simulations  show  that  the  equilibrium  solution  is  oscillatory  and 
undamped.  Figure  25  shows  the  oscillation  of  the  cathode  charge  flux;  Figure  26  is  a  plot  of 
the  evolution  of  the  potential  minimum.  Both  results  show  that  the  particles  contributing  to  the 
cathode  charge  flux  are  returned  at  a  point  which  varies  in  time. 

Figures  27  and  28  are  similar  diagnostics,  but  for  7V0,  T  small.  We  see  that,  for 
vint h/  vinbim  small  enough,  we  still  have  an  oscillatory  (although  initially  damped)  solution.  The 
damping  and  its  duration  increase  with  the  vinth/ vinbim  ratio.  On  a  series  of  runs  with  r—0 
and  T  small,  it  has  been  found  that  the  frequency  of  the  oscillation,  whether  initially  damped 
or  not,  varies  linearly  with  t»Ax+_J,  i.e.  with  the  plasma  frequency  at  the  point  (the  position 
of  which  is  time-dependent)  where  particles  are  returned  to  the  cathode. 

Now,  as  we  further  increase  the  temperature,  damping  becomes  stronger,  until  there  is 
virtually  no  more  oscillation  at  all.  Note  that  a  damped  oscillation  can  be  obtained  by  giving 
Wrtfli— 0.0,  but  nenter*  1:  even  though  all  particles  should  ideally  be  injected  at  the  same  velo¬ 
city,  since  vmrh— 0.0,  they  all  have  different  real  injection  velocities  (i.e.  their  velocities  depend 
on  their  positions  when  they  are  counted  inside  the  diode  for  the  first  time;  but  since 
ntnter*  1,  all  entering  particles  at  any  given  time  step  are  at  different  positions  (according  to 
Section  3.2.2),  i.e.  they  behave  as  if  they  had  different  real  injection  velocities.  Conclusion:  one 
should  only  use  nenter~l  when  one  wants  to  simulate  a  really  cold  case! 

Finally,  Figures  29  and  30  show  the  boundary  charge  fluxes  and  the  motion  of  4>mm  for  a 
warm  case  with  no  drift  in  the  reservoir.  One  has  then  a  purely  time-independent  solution  (or, 
to  be  more  accurate,  what  is  left  of  time-dependence  is  due  to  noise  and  to  the  numerical 
effects  described  in  Sections  4.1.2  and  4.1.3). 


10  PUSMA  SIMULATION  01/20/82  22:32:  1  1 
IOX  84 4  FIGURE  20 


HIPRI0 

m 

0.87, 

IE 

■1 

0, 

IFVX 

• 

0, 

IPACK  - 

B. 

1  PH  1 

m 

0, 

IRHO 

■ 

0, 

IRHOS 

m 

0, 

ivxvr  - 

0, 

IXVX 

m 

0. 

J ANODE 

■ 

•FALSE, 

JCATOO 

m 

.FALSE, 

JEM  1 TR  - 

.FALSE 

JT0TAL 

m 

.FALSE, 

NPRIO 

■ 

0. 

NRANK 

m 

3, 

MPL0T 

|EN0 

m 

0, 

A  1 

m 

0  D, 

A2 

m 

o.o, 

OT 

m 

0.01  , 

EPSI 

1.0, 

IW 

m 

2, 

L 

m 

3.141593 

,NG 

m 

312  , 

NSP 

1. 

NT 

m 

500, 

RANOVX 

m 

0.0, 

VEC 

m 

•TRUE, 

PH  1 0 

■40,2.5, 

0,  11. 

.0, 

5.0.0, 

$END 

MODE 

m 

1  , 

N 

■ 

0. 

NENTER 

m 

512, 

NMAX 

131072 

CM 

-to. 

THETAV 

m 

o.o. 

THETAX 

m 

0.0, 

VINBIM  - 

0.0, 

V  INMAX 

m 

1  0E+1O, 

V 1 NM 1 N 

m 

o.o. 

VINTH 

u 

1.0. 

VT1  ■ 

o.o, 

YT2 

m 

0.0, 

VO 

m 

o.o. 

VI 

m 

o.o, 

WC 

o.o, 

WP 

m 

0.1, 

XI 

m 

0.0, 

{END 

PLASMA 

SHEATH  SIMULATION 

IN  ID 

OIM*UBSNvaBBNVWSaN*ltBaNV«  as 


TOTAL  ENERGY,  TIME-  0.0000  TO  3.0000 


Figure  23  (Section  4.2.2):  Total  energy  during  settlement  of  two  successive  equilibria 
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Figure  24  (Section  4.2.2):  frmm  and  x*  during  settlement  of  equilibrium. 
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Figure  26  (Section  4.2.3):  Evolution  of  in  a  cold,  space-charge-limited  diode. 
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Figure  27  (Section  4.2.3):  Boundary  charge  fluxes  in  a  space-charge-limited  diode  with 

r*  0, 7*0. 


10  PLASMA  SIMULATION  08/23/12  21; 30:50 

SOX  94* 

F IGURE 

28 

HIPRI0 

m 

0.87, 

IE 

■ 

o. 

irvx 

a 

0, 

IPACK 

a 

8. 

1  PM  1 

m 

0, 

IRHO 

* 

0, 

IRHOS 

a 

o, 

IYXVY 

a 

0, 

IXVX 

m 

0, 

J ANODE 

■ 

.FALSE, 

JCATOO 

a 

.FALSE, 

JEM  1 TR 

a 

.false 

i TOTAL 

m 

■FALSE, 

RPR  1 0 

a 

0, 

DRANK 

a 

9. 

MPLOT 

m 

0, 

|EN0 

A 1 

m 

00, 

A2 

a 

0.0, 

OT 

a 

0.005, 

EPSI 

a 

1.0, 

IW 

m 

2 1 

L 

a 

3. 141 StS 

,NG 

a 

512  . 

NSP 

a 

1, 

Ml 

m 

*000, 

R ANDY  X 

a 

o.o, 

YEC 

a 

.TRUE, 

PH  10 

m 

0.0,00, 

0, 

(END 

MODE 

m 

r , 

N 

a 

0. 

RENTER 

a 

1024, 

NMAX 

a 

131072 

QM 

—  1  0, 

THETAV 

a 

o.o, 

1HETAX 

a 

o.o, 

VIKBIM 

a 

10. 0, 

Y INMAX 

■ 

1  0E+1Q, 

V 1 NM 1 N 

a 

0.0, 

VINTH 

a 

0.5, 

yrt 

a 

0.0, 

m 

• 

0.0, 

yo 

a 

0.0, 

VI 

a 

o.o, 

« 

a 

o.o, 

WP 

m 

1.0. 

XI 

a 

0.0, 

)EN0 

PLASMA 

SHEATH  SIMULATION 

i 

1  10 

ml 


MnSilHiGISi 


illSSI3i;=! 

!' 11 11  IIP  Si' 

I  :  Ml  ■'tili 


dill 


hill! 


PoaHfon  x 


3  3  3  3 

Po*W*n 


PNtRIIAl  MIRIHIR,  TIRE-  4. ION  TO  2.4010 


fOTtNT  1*1  MIMIRIR,  TIRE-  2. SOM  tO  J.< 


Figure  28  (Section  4.2.3):  Evolution  of  4>„u0  in  a  space-charge-limited  diode  with  r*0,  T* 0. 
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Figure  29  (Section  4.2.3):  Boundary  charge  fluxes  in  a  warm,  space-charge-limited  diode. 
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Figure  30  (Section  4.2.3):  Evolution  of  in  a  warm,  space-charge-limited  diode. 
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4.2.4.  Effect  of  on  the  Velocity  Distribution 

When  $o  increases,  less  and  less  particles  will  be  reflected  back  to  the  cathode,  and  the 
average  velocity  will  increase.  More  high-energy  particles  will  be  present  in  the  system  at  each 
instant. 

Figure  31  shows  the  Maxwellian  distribution  of  velocities  for  a  totally-reflecting  4>o 

In  figure  32,  <t>0  is  much  larger  than  in  the  previous  case;  this  can  be  deduced  from  the 
presence  of  a  hot  tail  in  the  distribution  of  velocities. 

This  confirms  the  obvious  fact  that  particles  receive  energy  during  their  transit  time  in  the 
diode,  provided  the  balance  current  is  nonzero  ^ 0). 


4.2.5.  /—/(♦  o)  Characteristic 

The  S-shaped  curve  giving  /  as  a  function  of  <t>0  is  well-known:  for  4>0  small  enough, 
4<w*—0;  then,  as  d>0  increases;  iaMdt  becomes  nonzero,  and  increases  with  <f>o,  until  it  reaches  a 
saturation  value  above  which  w. 

Figure  33  shows  iano<*  (<ho)  and  (4>o).  The  high  noise  levels  and  the  almost  straight 
shapes  of  the  currents  are  due  to  the  fact  that,  for  graphics  purpose,  we  let  d>0(r)  increase  fairly 
rapidly  in  time,  whereas  the  theoretical  S-curve  is  obtained  for  d>o^d>o(/)-  Also,  at„  is  time- 
varying,  since  the  average  density  within  the  sheath  varies;  this  means  that  on  this  plot  <t>0  does 
not  increase  linearly  with  t. 

Other  simulation  runs  made  on  larger  time-scales  have  given  much  'cleaner*  results, 
agreeing  quite  well  with  theory. 
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Figure  32  (Section  4.2.4):  /( vx)  in  the  diode  for  $o  reflecting  part  of  the  particles. 


10  PLASMA  SIMULATION  OB/23/B2  21:21:40 
BOX  844  FIGURE  33 


HIPRIO 

m 

0.87, 

IE 

■ 

0, 

irvx 

■ 

0, 

1  PACK  - 

B, 

1  PH  1 

m 

0, 

IRHO 

■ 

0, 

IRHOS 

■ 

0, 

IVXVY  - 

0, 

IXVX 

m 

0, 

J ANODE 

■ 

•FALSE, 

JCATOO 

■ 

.FALSE, 

JEM  1  TR  - 

.FALSE 

J TOTAL 

MPLOT 

(END 

m 

m 

.FALSE, 

0, 

NPRIO 

■ 

0. 

NRANK 

■ 

3. 

A 1 

m 

o.o. 

A2 

■ 

0.0, 

OT 

• 

0.01 , 

EPSI 

1.0, 

IW 

m 

2, 

L 

■ 

3.14-15*3 

,  HO 

• 

312, 

NSP 

1, 

MT 

m 

19000, 

RANOVX 

• 

0.0, 

VEC 

■ 

.TRUE, 

PH  1 0 

■4.0,90, 

0,  11 

0 

190.0,1, 

$END 

MODE 

m 

1  , 

N 

• 

0. 

RENTER 

■ 

512, 

NMAX  - 

131072 

QM 

■1.0. 

THCTAV 

■ 

0.0. 

THETAX 

■ 

0.0, 

VtNBIM  - 

0.0, 

V  INMAX 

m 

1  .  DE+10, 

V 1 NM 1 N 

■ 

o.D, 

VINTH 

m 

1.0, 

VT1  - 

o.o, 

VT2 

m 

00. 

VO 

■ 

0.0, 

VI 

m 

0-0. 

fC 

o.o, 

WP 

m 

01. 

XI 

■ 

0.0, 

(END 

PLASMA  SHEATH  SIMULATION  IN  10 


O-FLUXES  AT  BOUNOAN IES  (SPECIES  1),  TIME-  0.0000  TO  190.0000 

Figure  33  (Section  4.2.5):  *-/(♦  o). 


-45- 


iso 


Time*  Dependent  Chi  Id- Langmuir  Diode  Simulation 


-46- 


5.  CONCLUSION 

We  believe  that  the  goal  of  this  project  has  been  reached,  since  we  have  now  a  tool  which 
allows  us  to  better  understand  a  wide  range  of  physical  phenomena  in  several  problems  related 
to  the  plasma  sheath. 

In  the  future,  we  expect  other  people  and  ourselves  to  use  our  code  for  the  study  of  more 
elaborated  problems  involving  several  species,  a  magnetic  field,  possible  secondary  emission  at 
the  target,  etc. 

We  hope  that  the  basic  information  contained  in  this  report  will  be  beneficial  in  helping 
those  who  wish  to  understand  our  code  and/or  to  use  it. 
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